Formation R initiation

Qu’est ce que R ?

R (https://www.r-project.org/) est un langage de programmation dédié à l’analyse statistique. Libre et gratuit, disponible sur toutes les plateformes (Mac OS, Windows, Linux), il s’est largement imposé ces dernières années dans le domaine des sciences humaines. R dispose par défaut d’instructions de base permettant de réaliser les opérations les plus courantes. De nouvelles fonctionnalités, regroupées en packages, peuvent être facilement ajoutées pour venir enrichir les possibilités du langage. Le point fort de R par rapport à des logiciels à interfaces graphiques, est qu’il permet très facilement l’automatisation et la reproductibilité des traitements.

Environnement de travail

Rstudio (https://www.rstudio.com/) est un environnement de travail permettant d’écrire du code R. Il ne s’agit pas d’un environnement clic bouton mais d’une véritable interface de développement (IDE) destiné à écrire du code. Il s’agit en fait d’une couche logiciel permettant d’écrire, d’exécuter et de visualiser les résultats d’un programme écrit dans le langage R. RStudio n’est pas indispensable pour utiliser R. Néanmoins, il est très pratique et très bien fait. [En savoir plus].

Le menu Session > Set Working Directory permet de sélectionner le répertoire par défaut contenant vos scripts et éventuels fichiers de données. Vous pouvez également le modifier via une ligne de code :

setwd('C:/Users/Stephanie/Desktop/R_initiation')

Dans Rstudio pour compiler une ligne, vous pouvez placer le curseur dessus et cliquer sur Run, ou bien utiliser le raccourci clavier ctrl+entrée

Opérations de base

Sous R, les éléments de base sont des objets : vecteurs, matrices, listes, table … Ces objets peuvent contenir des éléments de type numérique, booléen (logical : TRUE, FALSE), chaîne de caractère (string), facteur (factor, pour les variables qualitatives prenant un nombre déterminé de modalité). La différence entre la table de données (data.frame) et la matrice est tient notamment dans le fait que la première peut contenir des éléments mixtes.

Objets et opérations simples

On commence par créer des éléments pour illustrer les opérations de base en R. Le signe <- permet d’assigner un contenu à une variable dont on indique le nom à gauche.

element_a <- 2
element_a

## [1] 2

element_b <- 3
element_b

## [1] 3

element_c <- element_a/2 + 2*element_b # les opérations -, *, /, sqrt pour racine, ^ pour puissance, log, exp, sont également possibles !
element_c

## [1] 7

On peut travailler très naturellement sur des vecteurs

vecteur_a <- c(2,3,5) 
vecteur_a

## [1] 2 3 5

On peut récupérer un élément d’un vecteur avec l’opérateur [

vecteur_a

## [1] 2 3 5

vecteur_a[1] #on commence à numéroter à partir de 1 (et non de zéro)

## [1] 2

vecteur_a[c(1,3)]

## [1] 2 5

Dans le cas où le vecteur correspond à une séquence de nombres, on peut utiliser une syntaxe particulière

vecteur_a <- c(2,3,4) 
vecteur_b <- 2:4
vecteur_c <- seq(2,4,1) # les arguments de la fonction seq (séquence) correspondent à from, to, by (à partir de 1, jusqu'à 10, avec des écarts de 2)

vecteur_a == vecteur_b # compare les éléments un à un

## [1] TRUE TRUE TRUE

variable_logique_a <- vecteur_a == vecteur_b

Les opérations sur les vecteurs sont très proches des opérations sur les éléments. Une opération entre un vecteur a et un vecteur b revient à réaliser l’opération entre les couples d’éléments issus des deux vecteurs et situés à la même place (il faut que les vecteurs admettent la même taille pour que ça ait un sens)

vecteur_c <- 2*vecteur_c # on écrase la valeur du vecteur_c
vecteur_c + vecteur_a/2 # exemple d'opération entre deux vecteurs

## [1]  5.0  7.5 10.0

Pour concaténer deux vecteurs, on fera tout simplement

vecteur_a

## [1] 2 3 4

vecteur_b

## [1] 2 3 4

vecteur_d <- c(vecteur_a,vecteur_b)
vecteur_d

## [1] 2 3 4 2 3 4

Le vecteur est caractérisé notamment par sa taille

length(vecteur_a)

## [1] 3

Il peut-être trié :

vecteur_d <- c(3,2,1)
order(vecteur_d)

## [1] 3 2 1

sort(vecteur_d)

## [1] 1 2 3

Remarque, si on ne sait plus la signification des trois arguments de seq, on peut aller dans l’aide avec ?seq

liste_a <- list(2,3,5)
liste_a

## [[1]]
## [1] 2
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 5

liste_a[[2]]

## [1] 3

liste_b <- list(sexe = 2, age = 3, salaire = 5)
liste_b$salaire # si on donne des noms aux éléments de la liste, on peut les récupérer via l'opérateur $

## [1] 5

Pour concaténer deux listes

liste_a

## [[1]]
## [1] 2
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 5

append(liste_a, element_a) # ou liste_b à la place de element_a pour concaténer deux listes

## [[1]]
## [1] 2
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 5
## 
## [[4]]
## [1] 2

liste_a

## [[1]]
## [1] 2
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 5

La fonction length convient aussi pour les listes

length(liste_a)

## [1] 3

on peut passer de la liste au vecteur avec la fonction unlist

unlist(liste_a)

## [1] 2 3 5

Pour les matrices :

matrice_a <- matrix(1:15,ncol=5) # exemple d'une matrice remplie des chiffres consécutifs de 1 à 15 rangés sur 5 colonnes 
matrice_a

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    4    7   10   13
## [2,]    2    5    8   11   14
## [3,]    3    6    9   12   15

head(matrice_a) # la fonction head permet de visualiser un extrait des données, ici elles sont petites donc c'est la totalité

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    4    7   10   13
## [2,]    2    5    8   11   14
## [3,]    3    6    9   12   15

matrice_a[1,2] # pour récupérer l'élément de la première ligne et de la deuxième colonne 

## [1] 4

Par défaut, les opérations mathématiques simples se font élément par élément

matrice_b <- 2*matrice_a
matrice_a

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    4    7   10   13
## [2,]    2    5    8   11   14
## [3,]    3    6    9   12   15

matrice_b

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    2    8   14   20   26
## [2,]    4   10   16   22   28
## [3,]    6   12   18   24   30

matrice_b - matrice_a # les opérations se font éléments par élément

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    4    7   10   13
## [2,]    2    5    8   11   14
## [3,]    3    6    9   12   15

t(matrice_b)%*% matrice_a # pour faire des vrais produits matriciels, on utilise l'opérateur %*%, ici t() indique que l'on prend également la transposée

##      [,1] [,2] [,3] [,4] [,5]
## [1,]   28   64  100  136  172
## [2,]   64  154  244  334  424
## [3,]  100  244  388  532  676
## [4,]  136  334  532  730  928
## [5,]  172  424  676  928 1180

Pour concaténer deux matrices, on peut soit les mettre côte à côte (cbind), soit les empiler (rbind). Ces fonctions fonctionnent aussi pour les data.frames.

matrice_c <- cbind(matrice_a, matrice_b)
matrice_c

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    1    4    7   10   13    2    8   14   20    26
## [2,]    2    5    8   11   14    4   10   16   22    28
## [3,]    3    6    9   12   15    6   12   18   24    30

dim(matrice_c)

## [1]  3 10

matrice_d <- rbind(matrice_a, matrice_b)
matrice_d

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    4    7   10   13
## [2,]    2    5    8   11   14
## [3,]    3    6    9   12   15
## [4,]    2    8   14   20   26
## [5,]    4   10   16   22   28
## [6,]    6   12   18   24   30

dim(matrice_d)

## [1] 6 5

Les dimensions de la matrice peuvent être obtenues de la manière suivante :

nrow(matrice_a)

## [1] 3

ncol(matrice_a)

## [1] 5

dim(matrice_a)

## [1] 3 5

Les tables de données ou data.frame est sans doute le format qu’on utilisera le plus dans la suite.

df_a <- as.data.frame(matrice_a) # on commence par transformer la matrice en data.frame, l'opération symétrique as.matrix() est possible
head(df_a) 

names(df_a) # ici le data.frame provient de la conversion d'une matrice, il n'y a donc pas de noms de colonnes à part ceux qui sont donnés par défaut V1, V2...

## [1] "V1" "V2" "V3" "V4" "V5"

names(df_a) <- c('a','b','c','d','e') # on peut changer le nom des colonnes simplement

Pour récupérer de l’information dans la table de données, on peut procéder de différentes manières

df_a$b # avec l'opérateur $, on peut récupérer directement la colonne b

## [1] 4 5 6

df_a[['b']] # revient au même

## [1] 4 5 6

df_a[[2]] # correspond également car il s'agit de la 2e colonne

## [1] 4 5 6

df_a$e <- 1 # permet de construire une nouvelle colonne e remplie de 1 
head(df_a)

df_a$f <- df_a$a + df_a$b # permet de construire une nouvelle colonne qui serait la somme des deux premières, les opération -, *, /, sqrt pour racine carrée, log, exp, ^, sont également possibles ! 
head(df_a)

df_a$f <- sqrt(df_a$a) + df_a$b^2 # permet de construire une nouvelle colonne qui serait la somme des deux premières, les opération -, *, /, sqrt pour racine 

Vérifions le type de ces objets avec la fonction class

class(element_a)

## [1] "numeric"

class(vecteur_a)

## [1] "numeric"

class(variable_logique_a)

## [1] "logical"

class(liste_a)

## [1] "list"

class(matrice_a)

## [1] "matrix"

class(df_a)

## [1] "data.frame"

Nous avons créé un certain nombre de variables qui sont disponibles dans votre environnement (généralement en haut à droite par défaut) On peut aussi utiliser la fonction ls pour voir la liste des objets dont on dispose, et rm pour en supprimer.

ls()

##  [1] "df_a"               "element_a"          "element_b"         
##  [4] "element_c"          "liste_a"            "liste_b"           
##  [7] "matrice_a"          "matrice_b"          "matrice_c"         
## [10] "matrice_d"          "variable_logique_a" "vecteur_a"         
## [13] "vecteur_b"          "vecteur_c"          "vecteur_d"

rm(element_a)
ls()

##  [1] "df_a"               "element_b"          "element_c"         
##  [4] "liste_a"            "liste_b"            "matrice_a"         
##  [7] "matrice_b"          "matrice_c"          "matrice_d"         
## [10] "variable_logique_a" "vecteur_a"          "vecteur_b"         
## [13] "vecteur_c"          "vecteur_d"

Boucle, condition

Reprenons le cas où on a un data.frame simple. On veut calculer la moyenne, la somme, le maximum et le minimum de chaque colonne.
Une manière de le faire serait de faire une boucle sur le nombre de colonnes et d’appliquer la fonction mean, sum, max, min, à chaque colonne tour à tour.

nb_col <- ncol(df_a) 
mean_col <- NULL # on commence par créer une variable vide, dans laquelle on va ajouter itérativement les moyenne de chaque colonne

for (i in 1:nb_col){
  print(i)
  mean_temp <- mean(df_a[[i]])
  print(mean_temp)
  mean_col <- c(mean_col, mean_temp)
  print(mean_col)
  print('---------')
}

## [1] 1
## [1] 2
## [1] 2
## [1] "---------"
## [1] 2
## [1] 5
## [1] 2 5
## [1] "---------"
## [1] 3
## [1] 8
## [1] 2 5 8
## [1] "---------"
## [1] 4
## [1] 11
## [1]  2  5  8 11
## [1] "---------"
## [1] 5
## [1] 1
## [1]  2  5  8 11  1
## [1] "---------"
## [1] 6
## [1] 27.04875
## [1]  2.00000  5.00000  8.00000 11.00000  1.00000 27.04875
## [1] "---------"

En réalité, les boucles sont à éviter en R, elles ne sont pas efficaces, on leur préférera la fonction lapply qui distribue la fonction souhaitée sur chaque élément d’une liste. Un data.frame peut-être vu comme une liste de vecteurs.

mean_col <- lapply(df_a, mean) # ici la fonction est particulièrement simple car mean est déjà une fonction R
mean_col <- sapply(df_a, mean) # sapply procède de même mais renvoie un vecteur
mean_col

##        a        b        c        d        e        f 
##  2.00000  5.00000  8.00000 11.00000  1.00000 27.04875

Admettons qu’on cherche à appliquer une fonction personnalisée

mean_personnalisee <- function(vect){
  return(sum(vect)/length(vect))
}

mean_col <- sapply(df_a, function(vect) mean_personnalisee(vect))
mean_col

##        a        b        c        d        e        f 
##  2.00000  5.00000  8.00000 11.00000  1.00000 27.04875

Appliquée sur un vecteur, lapply le convertit en liste, on peut repasser à un format vecteur avec la fonction unlist

fonction_personnalisee <- function(x){
  return(x*log(x))
}

vecteur_a_transforme <- lapply(vecteur_a, function(x) fonction_personnalisee(x))
class(vecteur_a_transforme)

## [1] "list"

vecteur_a_transforme

## [[1]]
## [1] 1.386294
## 
## [[2]]
## [1] 3.295837
## 
## [[3]]
## [1] 5.545177

unlist(vecteur_a_transforme)

## [1] 1.386294 3.295837 5.545177

vecteur_a*log(vecteur_a) # aurait marché aussi !

## [1] 1.386294 3.295837 5.545177

apply s’applique aussi à des matrices, on indique alors si la fonction à distribuer doit être distribuée en ligne ou en colonne

matrice_a

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    4    7   10   13
## [2,]    2    5    8   11   14
## [3,]    3    6    9   12   15

apply(matrice_a,1,sum) # ici on applique la fonction somme, autrement dit on somme les éléments par ligne

## [1] 35 40 45

apply(matrice_a,2,sum) # par colonne

## [1]  6 15 24 33 42

La syntaxe de la condition est très proche de celle de la boucle, sauf que l’on remplace for par if

fonction_personnalisee <- function(x){
  if (x>0){res <- x*log(x)}
  else {res <- 0}
  return(res)
}

vecteur_d <- c(0,-1,2)
sapply(vecteur_d, function(x) fonction_personnalisee(x))

## [1] 0.000000 0.000000 1.386294

abs(vecteur_d) # rmq : abs permet de passer un élément en valeur absolue (s'applique également aux vecteurs sans recours à lapply)

## [1] 0 1 2

Ici on a inclus la condition dans une fonction, mais ce n’est pas indispensable ! On peut vouloir aussi tester la différence != plutôt que l’égalité, ‘<’ ou ‘>’, et si on a plusieurs conditions, on mettra chacune entre parenthèse et on utilisera | pour dire “ou” et & pour dire “et” (par exemple (b-3==a) & (b>=a)).

Parfois, on n’a pas besoin de passer par if pour appliquer une condition, par exemple si on veut récupérer seulement les éléments de vecteur_d supérieurs ou égaux à 0. Cela revient à dire qu’on teste la condition sur chaque élément et qu’on conserve vecteur_d[i] pour i tel que vecteur_d[i]>=0

vecteur_d

## [1]  0 -1  2

vecteur_d[vecteur_d>=0]

## [1] 0 2

Et ça fonctionne aussi avec les data.frames

df_a[df_a$a>2] # on ne garde que les colonnes qui vérifient cette condition

Le type chaîne de caractère

Il arrive que les objets que l’on manipule ne soit pas numériques, on parle alors de character ou de factor. C’est le cas par exemple des noms de colonnes de df_a

class(names(df_a))

## [1] "character"

On peut aussi réaliser des opérations, par exemple de concaténation. Imaginons que l’on souhaite préciser les noms de colonnes en indiquant qu’il s’agit de département.

paste('departement','test',sep='_')

## [1] "departement_test"

paste(rep('departement',5), names(df_a),sep='_') 

## [1] "departement_a" "departement_b" "departement_c" "departement_d"
## [5] "departement_e" "departement_f"

names(df_a) <- paste(rep('departement',5), names(df_a),sep='_') # revient à faire apply de paste sur tous les éléments du vecteur names(df_a)
rep('departement',5) # rep est la fonction qui permet de créer des vecteurs d'une taille donnée contenant toujours le même élément qui pourrait être numérique

## [1] "departement" "departement" "departement" "departement" "departement"

On peut récupérer une partie d’une chaine de caractère simplement

character_a <- 'departement_a' 
class(character_a)

## [1] "character"

substr(character_a,1,10)

## [1] "departemen"

strsplit(character_a,'_')

## [[1]]
## [1] "departement" "a"

strsplit(character_a,'_')[[1]][[2]]

## [1] "a"

La taille d’une chaine de caractères correspond au nombre de caractères

nchar(character_a)

## [1] 13

Les facteurs correspondent aux modalités d’une variable qualitative.

col_df_a <- names(df_a)
col_df_a_fact <- as.factor(col_df_a)
is.character(col_df_a_fact)

## [1] FALSE

class(col_df_a)

## [1] "character"

class(col_df_a_fact)

## [1] "factor"

On ne peut pas rajouter un élément qui n’est pas dans la liste des facteurs naïvement

levels(col_df_a_fact)

## [1] "departement_a" "departement_b" "departement_c" "departement_d"
## [5] "departement_e" "departement_f"

col_df_a_fact <- c(col_df_a_fact, 'departement_z')
class(col_df_a_fact) # attention cette opération a converti col_df_a_fact en character

## [1] "character"

col_df_a_fact 

## [1] "1"             "2"             "3"             "4"            
## [5] "5"             "6"             "departement_z"

levels(col_df_a_fact) 

## NULL

col_df_a_fact <- as.factor(col_df_a)
as.character(col_df_a_fact) # parfois on préfère revenir aux chaines de caractères pour éviter ce type de problèmes 

## [1] "departement_a" "departement_b" "departement_c" "departement_d"
## [5] "departement_e" "departement_f"

Lorsque l’on charge des tables de données, il est possible que R ne reconnaisse pas une variable numérique et qu’elle soit codée en facteurs. Dans ce cas, il faut la convertir en chaines de caractères avant de la convertir en nombre, car les facteurs sont associés à des numéros qui correspondent aux modalités.

nombres <- c(4,2,3)
facteurs <- as.factor(nombres)
facteurs

## [1] 4 2 3
## Levels: 2 3 4

as.numeric(facteurs)

## [1] 3 1 2

as.numeric(as.character(facteurs))

## [1] 4 2 3

Générer des variables aléatoires

On a souvent besoin de générer des nombres aléatoires, par exemple pour tirer un échantillon dans les lignes d’une table, le plus simple est d’utiliser sample :

sample(1:100,5) # on tire 5 éléments sans remise entre 1 et 100, l'argument replace = TRUE permet de faire un tirage avec remise.

## [1] 87 17 16 31 62

Il est parfois intéressant de simuler des données. Générons 50 observations tirées d’une variable suivant une loi normale de moyenne 20 et d’écart-type 10.

set.seed(123) # permet de reproduire les mêmes résultats d'une fois sur l'autre en dépit de l'aléa présent
n <- 50
Y <- rnorm(n,20,5)

mean(Y) # moyenne

## [1] 20.17202

sd(Y) # écart-type

## [1] 4.62935

sd(Y)/sqrt(length(Y)) # écart-type de la moyenne

## [1] 0.6546889

summary(Y) # quartiles et moyenne

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.17   17.20   19.64   20.17   23.49   30.84

On peut représenter graphiquement cette variable

boxplot(Y) # diagramme boîte

hist(Y, probability=T, col="blue") # histogramme de la densité

lines(density(Y), col="red", lwd=2) # estimation par la méthode du noyau
# tracer la loi théorique
x <- 1:100
curve(dnorm(x,mean=20,sd=5),add=TRUE,col="green",lwd=2) # l'argument add permet de rajouter ce contour sur le graphique

On peut générer des observations tirées dans une loi uniforme sur [a,b]

X <- runif(n,50,100)
hist(X)

Vous pouvez augmenter la valeur de n pour voir ce que ça donne

Lecture et exploration de données

Pour cette entrée en matière de l’exploration d’un vrai fichier de donnée, nous nous sommes fortement inspirés de l’exercice suivant https://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-scenar-statlab.pdf

Une étude réalisée entre 1961 et 1973 dans la maternité d’un hôpital d’Oakland (Californie) avait pour but de rechercher si certaines caractéristiques des parents avaient une influence sur le développement de l’enfant. Parmi les variables collectées, 19 variables décrites dans le tableau ci-dessous ont été observées sur 115 familles ou unités statistiques. Ces variables décrivent des informations médicales et socio-économiques concernant le bébé et ses parents au moment de la naissance puis dix ans plus tard, permettant ainsi de se poser différentes questions de nature plutôt épidémiologique : - Influence ou non de la consommation de cigarettes sur le sexe de l’enfant, sur son poids, sur sa taille, - sur l’évolution du poids de la mère en 10 ans, - sur les liaisons entre les caractéristiques des parents (poids, taille, rhésus) et celles de leur enfant

Chargement et exploration sommaire

famil=read.csv2("statlab.csv") # lecture du fichier csv
head(famil) # aperçu du haut des données

dim(famil) # combien de colonnes et de ligne 

## [1] 115  19

summary(famil) # quelques statistiques classiques

##  ESx    ERh          ET0             EP0             ET10      
##  F:46   R-:22   Min.   :43.00   Min.   :1.770   Min.   :124.0  
##  M:69   R+:93   1st Qu.:51.00   1st Qu.:2.925   1st Qu.:132.5  
##                 Median :52.00   Median :3.360   Median :136.0  
##                 Mean   :51.93   Mean   :3.389   Mean   :136.0  
##                 3rd Qu.:53.50   3rd Qu.:3.760   3rd Qu.:139.8  
##                 Max.   :58.50   Max.   :6.350   Max.   :151.5  
##       EP10       MRh          MA0             MP0             MCig0   
##  Min.   :23.10   R-:19   Min.   :15.00   Min.   : 43.50   >10cig :39  
##  1st Qu.:27.45   R+:96   1st Qu.:24.00   1st Qu.: 51.50   0cig   :48  
##  Median :31.80           Median :27.00   Median : 56.50   1-10cig:28  
##  Mean   :32.53           Mean   :27.53   Mean   : 59.57               
##  3rd Qu.:36.05           3rd Qu.:32.00   3rd Qu.: 63.25               
##  Max.   :57.20           Max.   :43.00   Max.   :113.50               
##        MT             MP10            MCig10        PA0       
##  Min.   :148.5   Min.   : 38.50   >10cig :39   Min.   :20.00  
##  1st Qu.:160.2   1st Qu.: 55.25   0cig   :55   1st Qu.:26.00  
##  Median :164.0   Median : 61.00   1-10cig:21   Median :29.00  
##  Mean   :164.1   Mean   : 64.07                Mean   :30.57  
##  3rd Qu.:167.8   3rd Qu.: 69.25                3rd Qu.:35.00  
##  Max.   :178.0   Max.   :115.50                Max.   :50.00  
##      PCig0          PT             PP10             RF0       
##  >10cig :65   Min.   :158.0   Min.   : 49.50   Min.   : 14.0  
##  0cig   :36   1st Qu.:172.8   1st Qu.: 72.50   1st Qu.: 54.0  
##  1-10cig:14   Median :179.0   Median : 80.50   Median : 72.0  
##               Mean   :179.0   Mean   : 82.23   Mean   : 76.1  
##               3rd Qu.:185.5   3rd Qu.: 90.50   3rd Qu.: 93.0  
##               Max.   :198.0   Max.   :129.00   Max.   :250.0  
##       RF10      
##  Min.   : 21.0  
##  1st Qu.:110.0  
##  Median :150.0  
##  Mean   :156.5  
##  3rd Qu.:190.0  
##  Max.   :500.0

La fonction summary fournit déjà beaucoup d’information sur chaque variable. Mais on pourrait réappliquer les fonctions vues précedemment pour calculer ces indicateurs et/ou les représenter graphiquement (moyenne, écart-type, quantiles, diagramme boîte, histogramme). Cela permet notamment de repérer les valeurs atypiques, l’hétérogénéité des variances, les distributions dissymétriques…

sapply(famil, mean) # les moyennes de chaque variable 

##        ESx        ERh        ET0        EP0       ET10       EP10 
##         NA         NA  51.930435   3.388783 136.047826  32.526087 
##        MRh        MA0        MP0      MCig0         MT       MP10 
##         NA  27.530435  59.573913         NA 164.104348  64.065217 
##     MCig10        PA0      PCig0         PT       PP10        RF0 
##         NA  30.565217         NA 178.965217  82.234783  76.095652 
##       RF10 
## 156.495652

sapply(famil, sd) # les écarts type

##        ESx        ERh        ET0        EP0       ET10       EP10 
##  0.4920419  0.3950495  2.9793787  0.6767183  6.0244290  6.1476438 
##        MRh        MA0        MP0      MCig0         MT       MP10 
##  0.3730019  5.8524575 12.9283185  0.7605851  5.5025894 13.3862200 
##     MCig10        PA0      PCig0         PT       PP10        RF0 
##  0.7082385  6.4673681  0.7032669  7.7521140 13.6865476 36.1429845 
##       RF10 
## 72.5150165

boxplot(famil$ET0) # taille de l'enfant

boxplot(famil$EP0) # poids de l'enfant

hist(famil$ET0)

hist(famil$EP0)

Pour les variables qualitatives, on aura plutôt tendance à considérer les fréquences des modalités.

table(famil$ESx)

## 
##  F  M 
## 46 69

barplot(table(famil$ESx)) # sexe de l'enfant

barplot(table(famil$MCig0)) # consommation de cigarettes

pie(table(famil$MCig10)) # consommation de cigarettes dix ans après 

Pour analyser grossièrement les liaisons, on peut regarder un nuage de points

pairs(famil[,c(3:6,8,9,11,12,14,16:19)])

plot(EP10~PP10,data=famil) # poids de l'enfant à 10 ans, poids du père

plot(EP10~ET10,data=famil) # poids de l'enfant à 10 ans, taille à 10 ans

Pour analyser grossièrement les liaisons entre variables qualitatives, on aura recours aux tables de contingence

table(famil$ESx,famil$ERh) # sexe et rhésus

##    
##     R- R+
##   F  9 37
##   M 13 56

addmargins(table(famil$ESx,famil$ERh)) # avec les marges

##      
##        R-  R+ Sum
##   F     9  37  46
##   M    13  56  69
##   Sum  22  93 115

prop.table(table(famil$ESx,famil$ERh)) # fréquences relatives

##    
##             R-         R+
##   F 0.07826087 0.32173913
##   M 0.11304348 0.48695652

Une manière de représenter graphiquement ces tables est d’utiliser une moisaique

mosaicplot(table(famil$ESx,famil$ERh))

addmargins(table(famil$MCig0,famil$ESx)) # consommation de cigarette et sexe de l'enfant

##          
##             F   M Sum
##   >10cig   11  28  39
##   0cig     27  21  48
##   1-10cig   8  20  28
##   Sum      46  69 115

mosaicplot(table(famil$MCig0,famil$ESx))

Si l’on veut croiser variables qualitatives et quantitatives, on peut utiliser les boites

boxplot(EP0~ESx,data=famil) # poids de l'enfant vs sexe

boxplot(EP0~MCig0,data=famil) # poids de l'enfant vs consommation de cigarettes

Tests de liaison

Quelques tests, imaginons qu’on veuille tester l’indépendance de deux variables qualitatives. Le test du \(\Khi^2\) sera adapté à ce problème (dans le cas où les effectifs d’une modalité sont trop faibles, il faut regrouper les modalités).

chisq.test(table(famil$ESx,famil$ERh)) # Sexe de l'enfant vs Rhésus

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(famil$ESx, famil$ERh)
## X-squared = 5.6549e-31, df = 1, p-value = 1

chisq.test(table(famil$ESx,famil$MCig0)) # Sexe de l'enfant vs consommation de cigarettes

## 
##  Pearson's Chi-squared test
## 
## data:  table(famil$ESx, famil$MCig0)
## X-squared = 9.0657, df = 2, p-value = 0.01075

On s’intéresse à l’influence du sexe sur la taille à la naissance. Tester l’égalité des deux moyennes nécessite de vérifier préalablement plusieurs points : la normalité des distributions dans chaque classe à moins que l’échantillon soit considéré de taille suffisamment grande, le caractère indépendant ou appariés des échantillons, l’égalité ou non des variances à l’intérieur de chaque groupe. On dispose de deux échantillons indépendants : les garçons et les filles. Testons les autres hypothèses.

# Normalité des distributions (facultatif car n grand ici)
shapiro.test(famil[famil$ESx=="M","ET0"]) 

## 
##  Shapiro-Wilk normality test
## 
## data:  famil[famil$ESx == "M", "ET0"]
## W = 0.95676, p-value = 0.01779

shapiro.test(famil[famil$ESx=="F","ET0"])

## 
##  Shapiro-Wilk normality test
## 
## data:  famil[famil$ESx == "F", "ET0"]
## W = 0.95141, p-value = 0.05315

# égalité des variances (test de Fisher)
var.test(ET0~ESx,data=famil)

## 
##  F test to compare two variances
## 
## data:  ET0 by ESx
## F = 0.92761, num df = 45, denom df = 68, p-value = 0.7979
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5495201 1.6132113
## sample estimates:
## ratio of variances 
##          0.9276114

Le test de comparaison des moyennes à utiliser (Student vs. Welsh) dépend du résultat précédent concernant l’égalité des variances.

t.test(ET0~ESx,var.equal=F, data=famil) # si les variances sont différentes c'est un test de Welsh

## 
##  Welch Two Sample t-test
## 
## data:  ET0 by ESx
## t = 2.9017, df = 99.064, p-value = 0.004573
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.5006535 2.6660131
## sample estimates:
## mean in group F mean in group M 
##        52.88043        51.29710

t.test(ET0~ESx,var.equal=T, data=famil) # Dans le cas où elles sont considérées égales, c'est un test de Student.

## 
##  Two Sample t-test
## 
## data:  ET0 by ESx
## t = 2.8798, df = 113, p-value = 0.004761
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.4940799 2.6725868
## sample estimates:
## mean in group F mean in group M 
##        52.88043        51.29710

Dans le cas d’échantillons appariés, par exemple si on se propose d’étudier l’évolution du poids de la mère au moment de la naissance et dix ans après, on utilise l’option paired=TRUE

t.test(famil$MP0, famil$MP10,paired=TRUE)

## 
##  Paired t-test
## 
## data:  famil$MP0 and famil$MP10
## t = -7.458, df = 114, p-value = 1.864e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.684285 -3.298324
## sample estimates:
## mean of the differences 
##               -4.491304

Remarques : - si l’hypothèse de normalité des distributions n’est pas vérifiée et si l’échantillon est trop réduit, c’est un test non-paramétrique qu’il faut mettre en ouvre. Les tests non-paramétriques sont basés sur les rangs des observations et donc sur les comparaisons des médianes des échantillons (wilcoxson).

• Si l’on veut tester l’indépendance entre une variable qualitative et quantitative, l’ANOVA associée à un test de Fisher est sans doute le test le plus utilisé ; il revient au test de Student lorsque la variable qualitative n’a que deux modalités.

Régression simple et multiple

La régression simple permet de tester l’influence éventuelle d’une variable sur une autre et plus précisément, dans le cas de cet exemple, d’expliquer la taille de l’enfant à 10 ans. On estime le modèle avec la fonction lm

res1.reg=lm(ET10 ~ PT, data = famil)
names(res1.reg) # liste des résultats

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

Les graphiques suivants permettent de s’assurer de la validité du modèle, et notamment de statuer sur l’homoscédasticité des résidus, leur normalité, la bonne linéarité du modèle.

qqnorm(res1.reg$residuals) # normalité des résidus
qqline(res1.reg$residuals)

shapiro.test(res1.reg$residuals) 

## 
##  Shapiro-Wilk normality test
## 
## data:  res1.reg$residuals
## W = 0.98511, p-value = 0.2346

Les résidus sont “grands” si, une fois normalisés ou plutôt “studentisés”, ils sont de valeur absolue plus grande que 2.

res.student=rstudent(res1.reg)
ychap=res1.reg$fitted.values
plot(res.student~ychap,ylab="Résidus") 
abline(h=c(-2,0,2),lty=c(2,1,2)) # rajoute une ligne horizontale

Une observation est influente si elle a un grand résidu est est associée à une grande valeur sur la diagonale de la hat matrix. Cela correspond à une valeur élevée (plus grande que 1) de la distance de Cook.

cook=cooks.distance(res1.reg) # repérage de points influents
plot(cook~ychap,ylab="Distance de Cook")
abline(h=c(0,1),lty=c(1,2))

La significativité du modèle peut être analysée via la fonction summary :

summary(res1.reg)

## 
## Call:
## lm(formula = ET10 ~ PT, data = famil)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.2297  -3.8156  -0.1753   3.4956  12.5815 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 92.82881   12.44735   7.458 1.94e-11 ***
## PT           0.24149    0.06949   3.475 0.000725 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.751 on 113 degrees of freedom
## Multiple R-squared:  0.09657,    Adjusted R-squared:  0.08857 
## F-statistic: 12.08 on 1 and 113 DF,  p-value: 0.0007246

La régression linéaire simple conduit à un modèle très mal ajusté. Le modèle linéaire multiple ci-dessous, plus complexe, recherche un meilleur ajustement des données.

res2.reg=lm(ET10~ET0+EP0+MA0+MP0+MT+MP10+PA0+PT+PP10+RF0+RF10, data = famil)
plot(res2.reg) # diagnostics des résidus

summary(res2.reg) # résultats 

## 
## Call:
## lm(formula = ET10 ~ ET0 + EP0 + MA0 + MP0 + MT + MP10 + PA0 + 
##     PT + PP10 + RF0 + RF10, data = famil)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.753  -3.605  -0.583   3.575  13.055 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 19.395052  22.732841   0.853  0.39554   
## ET0          0.603298   0.293976   2.052  0.04269 * 
## EP0         -1.024211   1.330923  -0.770  0.44333   
## MA0         -0.070042   0.149033  -0.470  0.63937   
## MP0         -0.066802   0.082435  -0.810  0.41961   
## MT           0.346222   0.107974   3.207  0.00179 **
## MP10         0.078521   0.079884   0.983  0.32794   
## PA0          0.147236   0.130401   1.129  0.26148   
## PT           0.154598   0.089309   1.731  0.08644 . 
## PP10         0.035471   0.050199   0.707  0.48140   
## RF0         -0.007751   0.016938  -0.458  0.64819   
## RF10        -0.010480   0.008614  -1.217  0.22652   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.245 on 103 degrees of freedom
## Multiple R-squared:  0.3152, Adjusted R-squared:  0.242 
## F-statistic: 4.309 on 11 and 103 DF,  p-value: 2.792e-05

ACP

On peut aller plus loin dans l’exploration sur les données sans fixer une variable cible, en réalisant une analyse en composantes principales. On ne garde que les variables quantitatives.

Pour cela, on va faire appel à la librairie prcomp, ce n’est pas la seule qui permet de faire une ACP. Assez régulièrement sous R, on va avoir besoin de mobiliser des packages ou librairies. La liste complète des packages disponibles gratuitement est consultable sur le site du CRAN. L’installation d’un package supplémentaire peut se faire via le menu Packages>Installer le(s) package(s) et en choisissant un site miroir du CRAN. On peut également télécharger l’archive .zip correspondant au package et utiliser ensuite Packages/Installer depuis des fichiers zip

On peut sinon taper la commande :

#install.packages('prcomp') # décommenter pour installer

data=famil[,c(3:6,8,9,11,12,14,16:19)] # liste des varaibles quantitatives

noms=dimnames(data)[[2]];noms

##  [1] "ET0"  "EP0"  "ET10" "EP10" "MA0"  "MP0"  "MT"   "MP10" "PA0"  "PT"  
## [11] "PP10" "RF0"  "RF10"

names(data)

##  [1] "ET0"  "EP0"  "ET10" "EP10" "MA0"  "MP0"  "MT"   "MP10" "PA0"  "PT"  
## [11] "PP10" "RF0"  "RF10"

res.pca=prcomp(data,scale=T)
plot(res.pca) # décroissance des valeurs propres

summary(res.pca) # parts de variance expliquée

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6
## Standard deviation     1.8915 1.4792 1.3594 1.13565 1.04918 0.95572
## Proportion of Variance 0.2752 0.1683 0.1422 0.09921 0.08468 0.07026
## Cumulative Proportion  0.2752 0.4435 0.5857 0.68488 0.76956 0.83982
##                            PC7     PC8     PC9    PC10   PC11    PC12
## Standard deviation     0.82493 0.66919 0.56744 0.43219 0.4281 0.39240
## Proportion of Variance 0.05235 0.03445 0.02477 0.01437 0.0141 0.01184
## Cumulative Proportion  0.89217 0.92661 0.95138 0.96575 0.9798 0.99169
##                           PC13
## Standard deviation     0.32870
## Proportion of Variance 0.00831
## Cumulative Proportion  1.00000

biplot(res.pca) # biplot du premier plan principal

plot(res.pca$x,col=as.integer(famil$MCig0))
text(10*res.pca$rotation,noms,col="blue")
abline(h=0,v=0,lty=2)

Exploration avancée sur données des urgences

Cet exercice mobilise l’enquête nationale sur les structures des urgences hospitalières 2013. Un jour donné, le 11 juin 2013 de 8h à 8h le lendemain, un questionnaire papié renseigné en temps réel en parallèle de la prise en charge par tous les points d’accueils (736) des établissements de santé autorisés pour l’activité d’accueil et de traitement des urgences (634) et pour tous les patients (52018).

urgences <- read.csv2("enq_urgences_structure.csv", sep=',') # lecture du fichier csv
head(urgences) # aperçu du haut des données

dim(urgences) # combien de colonnes et de ligne 

## [1] 734 223

#names(urgences)
dim(urgences)

## [1] 734 223

Remarque, si le fichier était au format xls, on pourrait l’ouvrir, par exemple, avec le package xlsx

library("xlsx")
urgences2 <- read.xlsx("enq_urgences_structure.xls",1) # 1 correspond ici à l'index de la page 1 de l'excel
head(urgences2) # aperçu du haut des données

dim(urgences2)

## [1] 734 225

Pour chaque heure de pointage (8h, 12h, 18h, 22h, 8h) indiquée par les lettres A, B, C, D, E, on retient pour l’exercice les variables suivantes (cf questionnaire dans le dossier) :

Les variables de présence de patients aux urgences :

• 87 : Nombre de patients en cours d’évaluation, non admis (nombre de patients présents hors lit i.e. box, salle d’attente) pour lesquels une décision d’hospitalisation n’a pas été prise
• 88 : Nombre de patients admis et en attente d’hospitalisation sur un brancard ou un fauteuil roulant (« boarding »)
• 89 : Nombre de patients présents en lit (UHCD) qui ne sont pas en attente d’hospitalisation
• 90 : Nombre de patients présents en lit (UHCD) en attente d’hospitalisation

Les variables d’organisation du travail :

• 92 : Nombre de médecins urgentistes
• 93 : Nombre de médecins non urgentistes
• 94 : Nombre d’interne
• 95 : Nombre de médecins intérimaires et extérieurs
• 96 : Nombre de cadres de santé
• 97 : Nombre d’infirmiers diplômés d’État (IDE)
• 98 : Nombre d’aide soignants(AS)
• 99 : Nombre d’agents de service hospitalier (ASH)
• 100 : Nombre de secrétaires
• 101 : Nombre de brancardiers

Les variables liées à l’établissement :

A103 : Nombre de passages sur les 24 heures A4 : type d’accueil A2 : numéro finess de l’établissement

col_a_garder <- c(87:90,92:101)
nbcol <- length(col_a_garder)
col_a_garder <- unlist(lapply(c('A','B','C','D','E'), function(lettre) paste(rep('ID_',nbcol),rep(lettre,nbcol),col_a_garder,sep=''))) 
col_a_garder

##  [1] "ID_A87"  "ID_A88"  "ID_A89"  "ID_A90"  "ID_A92"  "ID_A93"  "ID_A94" 
##  [8] "ID_A95"  "ID_A96"  "ID_A97"  "ID_A98"  "ID_A99"  "ID_A100" "ID_A101"
## [15] "ID_B87"  "ID_B88"  "ID_B89"  "ID_B90"  "ID_B92"  "ID_B93"  "ID_B94" 
## [22] "ID_B95"  "ID_B96"  "ID_B97"  "ID_B98"  "ID_B99"  "ID_B100" "ID_B101"
## [29] "ID_C87"  "ID_C88"  "ID_C89"  "ID_C90"  "ID_C92"  "ID_C93"  "ID_C94" 
## [36] "ID_C95"  "ID_C96"  "ID_C97"  "ID_C98"  "ID_C99"  "ID_C100" "ID_C101"
## [43] "ID_D87"  "ID_D88"  "ID_D89"  "ID_D90"  "ID_D92"  "ID_D93"  "ID_D94" 
## [50] "ID_D95"  "ID_D96"  "ID_D97"  "ID_D98"  "ID_D99"  "ID_D100" "ID_D101"
## [57] "ID_E87"  "ID_E88"  "ID_E89"  "ID_E90"  "ID_E92"  "ID_E93"  "ID_E94" 
## [64] "ID_E95"  "ID_E96"  "ID_E97"  "ID_E98"  "ID_E99"  "ID_E100" "ID_E101"

col_a_garder <- c(col_a_garder, "ID_A103","ID_A4","ID_A2")
urgences <- urgences[,col_a_garder]
dim(urgences)

## [1] 734  73

Regardons d’abord la répartition des types d’établissement

table(urgences$ID_A4) # répartition des structures

## 
##                                                              Structures des urgences générales 
##                                                                                            536 
## Structures des urgences générales sur site ayant aussi une structure des urgences pédiatriques 
##                                                                                             94 
##                                                           Structures des urgences pédiatriques 
##                                                                                            104

On va regarder la distribution du nombre de passages total dans l’établissement, c’est un proxy de la taille de ce dernier. On va créer des modalités pour comparer des services d’urgences comparables. Pour les graphiques, on va désormais utiliser la librairie ggplot2, installez-là si elle n’est pas installée.

library(ggplot2)
ggplot(urgences, aes(x=ID_A103)) + geom_histogram(fill="blue", alpha=0.4) + ggtitle("distribution du nombre de passages") # alpha gère la transparence 

Comparons cette distribution par type de structure :

p <- ggplot(urgences, aes(x=ID_A103, fill = ID_A4)) + geom_histogram(alpha=0.4) + theme(legend.position="top") # en précisant la couleur par un nom de variables, on demande le tracé d'un histogramme pour chaque modalité 
p

Avec la fonction ggplotly de plotly on peut rendre le graphique interactif.

library(plotly)
ggplotly(p)

On va désormais transformer cette variable en variable qualitative pour rassembler les établissements par proxy de taille

valeurs_manquantes <- is.na(urgences$ID_A103) # on regarde si la variable admet des valeurs manquantes, 
head(valeurs_manquantes)

## [1] FALSE FALSE FALSE FALSE FALSE FALSE

(TRUE %in% valeurs_manquantes) # c'est le cas 

## [1] TRUE

sum(valeurs_manquantes) # nb de valeurs manquantes

## [1] 1

dim(urgences)

## [1] 734  73

urgences <- na.omit(urgences) # si on veut retirer toutes les lignes avec au moins une valeur manquante (pas toujours souhaitable)
dim(urgences)

## [1] 725  73

# on crée une nouvelle colonne
urgences$nb_passages_3 <- cut(urgences$ID_A103, breaks = c(0,40,80,max(urgences$ID_A103, na.rm=T))) # max, comme la plupart des fonctions renverra une erreur si appliqué sur un vecteur contenant des NA, on les enlève pour le calcul avec l'option na.rm 
table(urgences$nb_passages)

## 
##   (0,40]  (40,80] (80,263] 
##      187      309      229

On voudrait maintenant regarder le nombre de patients ou les effectifs par tranche horaire, on voit bien qu’il sera plus facile d’avoir une variable tranche horaire. Il suffit de récupérer les colonnes correspondant à chaque horaire, de les renommer, puis de les empiler en créant une colonne horaire

horaire_8 <- urgences[,c(paste('ID_','A',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
horaire_12 <- urgences[,c(paste('ID_','B',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
horaire_18 <- urgences[,c(paste('ID_','C',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
horaire_22 <- urgences[,c(paste('ID_','D',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
horaire_8_ <- urgences[,c(paste('ID_','E',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
head(horaire_8)

#on rajoute une colonne d'heure
horaire_8$heure <- '8h'
horaire_12$heure <- '12h'
horaire_18$heure <- '18h'
horaire_22$heure <- '22h'
horaire_8_$heure <- '8h_'

#pour pouvoir empiler les 5 tables il faut qu'elles aient les mêmes noms de colonnes 
names(horaire_8) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_12) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_18) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_22) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_8_) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')

urgences2 <- rbind(horaire_8, horaire_12, horaire_18, horaire_22, horaire_8_)
dim(urgences)

## [1] 725  74

dim(urgences2)

## [1] 3625   19

head(urgences2)

# on va remplacer les variables 'ID_A103', 'ID_A4' et 'ID_A2' par des noms de variables intelligibles
names(urgences2)[names(urgences2)=='ID_A103'] <- 'nb_passages'
names(urgences2)[names(urgences2)=='ID_A4'] <- 'type_structure'
names(urgences2)[names(urgences2)=='ID_A2'] <- 'index_structure'

On peut regarder le nombre de médecin urgentiste ID_92 en fonction de la tranche horaire et de la taille de la structure avec un graphique à boites.

sapply(urgences2, class) #attention certaines variables numériques ne sont pas considérées comme telles !

##           ID_87           ID_88           ID_89           ID_90 
##       "integer"       "integer"       "integer"       "integer" 
##           ID_92           ID_93           ID_94           ID_95 
##        "factor"        "factor"       "integer"        "factor" 
##           ID_96           ID_97           ID_98           ID_99 
##       "integer"        "factor"       "integer"       "integer" 
##          ID_100          ID_101     nb_passages  type_structure 
##       "integer"       "integer"       "integer"        "factor" 
##   nb_passages_3 index_structure           heure 
##        "factor"        "factor"     "character"

urgences2[,1:15] <- lapply(urgences2[,1:15],function(v) as.numeric(as.character(v)))

p<-ggplot(data = urgences2, aes(x=heure, y=ID_92, color=nb_passages_3)) + geom_boxplot() + ylab('Effectif de médecin urgentiste')
p

Si on veut plutôt regarder l’évolution du nombre moyen dans le cours de la journée, il va falloir d’abord calculer cette moyenne par groupe. Pour cela on va faire appel au package dplyr, vous devez l’installer si ce n’est pas déjà fait.

library(dplyr)

grouped <- group_by(urgences2, heure, nb_passages_3)
grouped_mean <- summarise(grouped, mean=mean(ID_92, na.rm = T))
head(grouped_mean)

# les deux actions peuvent être réalisées dans la foulée avec l'opérateur %>%
grouped_mean <- urgences2 %>% group_by(heure, nb_passages_3) %>% summarise(mean=mean(ID_92, na.rm = TRUE))
head(grouped_mean)

ggplot(data = grouped_mean, aes(x=heure, y=mean, group = nb_passages_3, color = nb_passages_3)) + geom_line() + geom_point()

Remarque : les fonctionnalités du packages dplyr sont bien plus riches, allez jeter un oeil à la cheatsheet data wrangling french fournie dans le dossier de la formation.

Maintenant on va regarder l’évolution de la distribution des effectifs de personnels et de patients dans le temps. Pour présenter les différents graphiques côte à côte, on va utiliser l’option facet de ggplot2 mais avant cela, il faut retravailler encore sur la forme des données. Commençons par regarder les distributions d’effectifs de professionnels, donc les variables de

urgences_ps <- urgences2[,c('index_structure','nb_passages_3','heure','ID_90','ID_92','ID_93','ID_94','ID_95','ID_96','ID_97','ID_98','ID_99','ID_100','ID_101')]
names(urgences_ps)[4:length(urgences_ps)] <-c('med urgentistes', 'med non urgentistes','internes','médecins int et ext', 'cadres','IDE','AS','ASH','secrétaires','brancardiers') 
library(reshape)
urgences_ps <- melt(urgences_ps, id=c("index_structure","nb_passages_3",'heure'))
head(urgences_ps)

names(urgences_ps)[4:5] <-c('PS','effectif')
head(urgences_ps)

On a donc créé une variable PS qui prend plusieurs modalités et on a pour chaque structure d’urgence identifiée par l’identifiant ID_A2, une ligne par profession avec les effectifs corresponds. Cette représentation rend la représentation graphique quasi immédiate avec ggplot2

urgences_ps$heure = factor(urgences_ps$heure , levels=c("8h", "12h", "18h", "22h", "8h_")) # on transforme la variable en facteur et on ordonne les heures
p <- ggplot(urgences_ps, aes(x=PS, y=effectif, fill=PS)) + geom_bar(stat = "identity") + facet_grid(as.factor(heure) ~ nb_passages_3) + theme(        axis.ticks.x=element_blank(),axis.text.x=element_blank())
p

# theme permet de retirer la légende sous l'axe des x, si vous retirez ce bout de code vous verrez la différence
# facet_grid permet de produire plusieurs graphique en fixant une variable ici le nombre de passage ou l'heure de pointage

Le même exercice peut être réalisé pour la distribution des types de patients.

Enfin, nous allons représenter la répartition des structures d’urgence en France. Pour cela nous allons utiliser la table sas ej.sas7bdat qui lie les codes FINESS des établissements avec leurs coordonnées géographiques. Cette base est stockée dans un format SAS, on va utiliser la librairie haven pour la lire.

library(haven)
ej <- read_sas(data_file = 'ej.sas7bdat')
head(ej)

Pour tracer sur une carte les structures d’urgence de notre table urgences, on va apparier notre table avec ej, sur la clé finess et ID_A2 respectivement. On va bien faire attention à garder tous les codes finess de la base urgences

urgences_loc <- merge(urgences[,c('ID_A2','ID_A4','ID_A103')],ej,by.x = 'ID_A2',by.y = 'finess', all.x = TRUE) #on garde aussi le type de structure et le nombre de passages
head(urgences_loc)

dim(urgences)

## [1] 725  74

dim(urgences_loc)

## [1] 725  15

dim(ej)

## [1] 55511    13

Le fait d’avoir indiqué l’option all.x = TRUE force à garder toutes les lignes de urgences même quand il n’existe pas de match possible. Ici certains codes finess ne semblent pas avoir de coordonnées géographiques disponibles. On ne va conserver que les structures qui ont des coordonnées

urgences_loc_sans_na <- na.omit(urgences_loc)
dim(urgences_loc_sans_na)

## [1] 329  15

Seulement la motié des structures sont géocodées.

On va regarder la répartition des structures d’urgence par type de structure (repéré par la variable ID_A4). La représentation de points sur une carte nécessite de maîtriser le concept de projection : https://medium.com/@_FrancoisM/introduction-%C3%A0-la-manipulation-de-donn%C3%A9es-cartographiques-23b4e38d8f0f.

library(leaflet)
library(rgdal)
urgences_loc_sans_na <- na.omit(urgences_loc)

# Ce bloc permet de convertir les coordonnées en Lambert 93 en WGS 84 plus traditionnel
coordinates(urgences_loc_sans_na) <- c("X", "Y")
proj4string(urgences_loc_sans_na) <- CRS("+init=epsg:2154") # WGS 84
CRS.new <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
urgences_loc_sans_na <- spTransform(urgences_loc_sans_na, CRS.new)
urgences_loc_sans_na$lon <- data.frame(coordinates(urgences_loc_sans_na))$X
urgences_loc_sans_na$lat <- data.frame(coordinates(urgences_loc_sans_na))$Y


factpal <- colorFactor(topo.colors(3), urgences_loc_sans_na$ID_A4) # palette de 3 couleurs

m <- leaflet(urgences_loc_sans_na) %>% addTiles() %>%
                         addCircles(lng = ~lon, lat = ~lat, weight = 1,
                                    radius = ~ID_A103*100, 
                                    popup = ~paste(rs, ", nb de passages : ", ID_A103, ", type de structure :", ID_A4),
                                    color = ~factpal(ID_A4), fillOpacity = 0.5) %>%
                         addLegend("bottomleft", title = "Type de structure d'urgences", pal = factpal, values = ~ID_A4, opacity = 0.7)
m